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ABSTRACT 

We study the Wouthuysen-Field couphng at early universe with numerical 
solutions of the integrodifferential equation describing the kinetics of photons 
undergoing resonant scattering. The numerical solver is developed based on the 
weighted essentially non-oscillatory (WENO) scheme for the Boltzmann-like in- 
tegrodifferential equation. This method has perfectly passed the tests of analytic 
solution and conservation property of the resonant scattering equation. We focus 
on the time evolution of the Wouthuysen-Field (W-F) coupling in relation to the 
21 cm emission and absorption at the epoch of reionization. We especially pay 
attention to the formation of the local Boltzmann distribution, ^-{'^~'^o)/kT ^ 
photon frequency spectrum around resonant frequency z/q within width z/;, i.e. 
I'^ — t'ol ^ ^i- We show that a local Boltzmann distribution will be formed if 
photons with frequency ~ vq have undergone a ten thousand or more times of 
scattering, which corresponds to the order of 10^ yrs for neutral hydrogen den- 
sity of the concordance ACDM model. The time evolution of the shape and 
width of the local Boltzmann distribution actually doesn't dependent on the de- 
tails of atomic recoil, photon sources, or initial conditions very much. However, 
the intensity of photon flux at the local Boltzmann distribution is substantially 
time-dependent. The time scale of approaching the saturated intensity can be as 
long as 10^-10^ yrs for typical parameters of the ACDM model. The intensity of 
the local Boltzmann distribution at time less than 10^ yrs is significantly lower 
than that of the saturation state. Therefore, it may not be always reasonable to 
assume that the deviation of the spin temperature of 21 cm energy states from 
cosmic background temperature is mainly due to the W-F coupling first stars 
or their emission/absorption regions evolved with a time scale equal to or less 
than Myrs. 
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- scattering 
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1. Introduction 

It is generally believed that the physical state of the universe at the epoch of reionization 
can be probed by detecting the redshifted 21 cm signals from the ionized and heated regions 
around the first generation of stars (e.g. Furlanetto, Oh, & Briggs 2006). The emission and 
absorption of 21 cm are caused by the deviation of the spin temperature Tg of neutral hydro- 
gen from the temperature of cosmic microwave background (CMB) Tcmb at the considered 
redshift z. Many calculations have been done on the 21 cm emission and absorption from 
ionized halos of the first stars (Chiizhoy et al. 2006; Cen 2006; Liu et al. 2007). A common 
assumption of these calculations is that the deviations of Tg from Tcmb are mainly due to 
the Wouthuysen-Field (W-F) coupling (Wouthuysen, 1952; Field, 1958, 1959). That is, the 
resonant scattering of Lya photons with neutral hydrogen atoms locks the color temperature 
Tc of the photon spectrum around the Lya frequency to be equal to the kinetic temperature 
of hydrogen gas T. Consequently, the spin degree of freedom is determined by the kinetic 
temperature of hydrogen gas T. 

The W-F coupling is from the kinetics of photons undergoing resonant scattering, which 
is described by a Boltzmann integrodifferential equation. All the above-mentioned calcula- 
tions are based on time-independent solution of the resonant scattering kinetic equation 
with Fokker-Planck approximation (Chen & Miralda-Escude, 2004; Hirata, 2006; Furlanetto 
& Pritchard 2006; Chuzhoy & Shapiro 2006). This is equal to assume that the time scale 
of the onset of the W-F coupling is less than all time scales related to the 21 cm emis- 
sion/absorption. However, even in the first paper of the W-F coupling, the problem of time 
scale has been addressed as follows; "One can infer from this fact that the photons (in a box), 
after an infinite number of scattering processes on gas atoms with kinetic temperature T, will 
obtain a statistical distribution over the spectrum proportional to the Planck-radiation spec- 
trum of temperature T. After a finite but large number of scattering processes, the Planck 
shape will be produced in a region around the initial frequency" (Wouthuysen, 1952). That 
is, the W-F coupfing is onset only "after a finite but large number of scattering". One 
cannot assume that the time-independent solution is available for the W-F coupfing if the 
time-scales of the evolution of the first stars and their emission/absorption regions are short. 
A study on the time evolution of radiation spectrum due to resonant scattering is necessary. 

This problem is especially important for the 21 cm signal from the first stars, as the life 
times of the first stars are short. The ionized and heated regions around the first stars are 
strongly time-dependent (Cen 2006; Liu et al. 2007). The 21 cm emission/absortion regions 
are located in a narrow shell just outside the ionized region. On the other hand, the speed of 
the ionization-front (I-front) is rather high, even comparable to the speed of light. The time 
scales of the formation and the evolution of the 21 cm regions can be estimated by ~ d/c. 
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d being the thickness of the shells of 21 cm emission and absorption. The time-independent 
solution would be proper only if Lya photons approach the time-independent state in a time 
shorter than that of the 21 cm region evolution. 

Very few works have been done on the time dependent behavior of the W-F coupling. 
There is the lack of a time dependent solution even for the Fokker-Planck approximation. 

The existed time-dependent solvers (e.g. Meiksin, 2006) cannot pass the tests of analytical 
solutions (Field, 1959). On the other hand, the WENO algorithm is found to be effective 
to solve the Boltzmann equation (Carrillo et al. 2003) and radiative transfer (Qiu et al. 
2006, 2007, 2008). In this paper, we will study the time-dependent behavior of the W-F 
coupling with the WENO method by numerically solving the integrodiffcrcntial equation. 
In this context, we also develop a numerical solver in accordance with the term of resonant 
scattering. 

The paper is organized as follows. Section 2 presents the basic equations of the res- 
onant scattering of photons. Section 3 very briefly mentions the numerical solver of the 
WENO scheme and its test with Field's analytic solutions, leaving the details of algorithmic 
issues to the Appendix. Section 4 presents the time-dependent W-F coupling with a static 
background. Section 5 shows the numerical results of the W-F coupling in an expanding 
background. Finally, conclusions are given in Section 6, in which the application to the 21 
cm problem is also addressed. 



2. Basic equations 

2.1. Radiative transfer equations with resonant scattering 

Since we focus on the time-scales of the W-F coupling, we consider a homogeneous and 
isotropically expanding infinite medium consisting of neutral hydrogen with temperature T. 
The kinetics of photons in the frequency space is described by the radiative transfer equation 
with resonant scattering (Hummer & Rybicki, 1992; Rybicki & Dell'antonio 1994) 

at vt ox 

-kc(f){x)J{x,t) + kc J n{ x,x')J{x',t)dx'-\-S{x,t) (1) 

where J is the flux in terms of the photon number in units s~^cm~^. H{t) = a{t)/a{t) is the 
Hubble parameter, a{t) being the cosmic factor; vt — i^hsT Irnf-I"^ is the thermal velocity 
of hydrogen atoms. The dimensionless frequency x is related to the frequency v and the 
resonant frequency vq\)Y x — {y — vq)/ /Si/d, and Ai/^) = uqVt/c is the Doppler broadening. 
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S{t,x) is the source of photons. The parameter k = x/Az/£), where x is the intensity of the 
resonant absorption given by x = Tie^riif 12/ mf.c, and rii being the number density of neutral 
hydrogen HI at ground state, /12 = 0.416 is the oscillator strength. The cross section of 
resonant scattering at the line center is 

ao = /i2(Az/,,)-\ (2) 

nieC 

In eq.([I]), (t){x) is the profile of the absorption line at the resonant frequency z/q. If the profile 
is dominated by Doppler broadening, we have 

= ^e--\ (3) 



The redistribution function 7l{x,x') of eq.([I]) gives the probability of a photon absorbed at 
frequency x', and isotropically re-emitted at frequency x. For coherent scattering, we have 
(Field 1958; Hummer 1962; Basko 1981) 

7^(x, x') = ^e2^^'+^'erfc[max(|x + b\, \x' + b\)], (4) 

where the parameter b = huo/rnvTC = 2.5 x 10~^(10^/T)^/^ is due to the recoil of atoms. It 
is in the range of 3 x 10~^ - 3 x 10"'^, if the temperature T is in the range of 1 K - 10^ K. 
The redistribution function is normalized as 

j n{x,x')dx = (j){x'). (5) 

Therefore, we have 

kc / (f){x)J{x,t)dx = kc / / Tl{x, x')J{x' ,t)dxdx' . (6) 



It means that the total number of photons absorbed given by the term kc(f){x)J{x, t) of eq.([T]) 
is equal to the total number of scattered photons. Therefore, with eq.([T]), the number of 
photons is conserved. 



2.2. Rescaling the equations 

We use the new time variable r defined as r = en 1 dot, which is in unit of the mean 
free flight time of photons at resonant frequency. For the concordance ACDM model, the 
number density of hydrog en atoms is nu = 1.88 x W-'^iVlbh'^ / 0.022) {1 + cm'^. Therefore, 
we have 

/ y N 1/2 . N 3 /o.o22\ 
^ = 0-054/hiM T7^ ^ 7VTT K ?/-^, (7) 



10 V Vi+'^y K^bh"^ J 
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where /hi = Jt-i/'^-h is the fraction of neutral hydrogen. 
We rescale the eq.([T]) by the following new variables 

J'{x,t) = [a{t)/a{0)]'^J{x,t), S'{x,t) = [a{t)/a{0)]^S{x,t). (8) 

Thus, eq.([T]) becomes 

dJ'jx, t) 

— q;^ — = -(j){x)J {x,t) 

r dv 

+ n{x,x')J'{x\T)dx' + + S'{x,t), (9) 
J ox 

where the parameter 7 is the so-called Sobolev parameter. 7 = {H/vrk) = {87iH/3Ai2^^ni) = 
{H m ei'Q / Tie^ riif 12) I where A is the wavelength for Lya transition. 7"^ measures the number 
of scattering during a Hubble time. It actually is the Gunn-Peterson optical depth given by 

Around the first stars, the number /hi is strongly dependent on time and position (Liu et 
al. 2007). It is as small as 10^^ within the ionized sphere, and as high as ~ 1 outside the 
ionized sphere. Therefore, the parameter 7 would be in the range from 1 to 10~^. 

The physical meaning of the terms on the right hand side of eq.([9]) is clear. The first 
term is the absorption at frequency x, the second term is the re-emission of photons with 
frequency x by scattering, and the third term describes the redshift of photons. The time 
scale of a photon moving Aa; in the frequency space is equal to 

Ar = 7~^Ax. (11) 

This actually is due to the Hubble expansion. 
Considering eq.([6]), eq.([9]) gives 

^ j J'dx = J S'dx. (12) 

This equation shows that the total number of photons J J'dx is dependent only on the 

sources, regardless of the parameter b of the resonant scattering. Since numerical errors 
accumulated over a long time evolution could be huge, Eq. f|T2|) is useful to check the reliability 
of a numerical code. We will use J for J' and S for S' in sections below. It will not cause 
confusion. 
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3. Numerical method 

We use the WENO scheme to solve the eq.Q. This algorithm has been given in Roy 
et al. (2009). Some of the algorithmic details are given in the Appendix. We present a test 
to show the good performance of our solver below. 

Figure 1 plots both the analytical (Field, 1959) and WENO numerical solutions of 
eq.([9]) with parameters 7 = and b = 0. It shows that the numerical solutions can follow the 
analytical solution in all the time t and the frequency x considered. This result is not trivial 
if compared with the results of other numerical solvers, such as Meiksin (2006), which shows 
a large deviation between the analytical and numerical solutions. Therefore, our scheme is 
more reliable. 



t= 10000 




-4 -2 2 4 

X 



Fig. 1. — Analytical (dashed) and WENO numerical (solid) solutions of eq.Q with 7 = 
and 6 = 0. The source is taken to be 5* = 0(x) and initial condition J(x, 0) = 0. 



4. Wouthuysen-Field coupling in a static background 

To study the effect of atomic recoil, we first solve the time evolution of J{x, r) in a 
static background, i.e. 7 = 0. A typical time-dependent result is shown in Figure 2, in 
which b = 0.03. The solution of Figure 2 is actually similar to that shown in Figure 1. 
Figure 1 shows that a fiat plateau around x = is to be formed at r > 100, while Figure 2 
shows that J{x, r) evolves into a Boltzmann distribution around x = as 

J{x, t) ~ J(0, r)e-2^^ = J(0, r)e-^('^-'^o)/'=^, |x| < xi. (13) 

This is the so-called "Planck shape in a region around the initial frequency" (Wouthuysen, 
1952). The expression eq. f|T3l) has also been found by Field (1959). We will call this feature 
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to be a local Boltzmann distribution. The width xi of the local Boltzmann distribution 
is numerically defined by the frequency range |x| < xi, in which the slope In J{x,T)/dx 
deviating from 2b is small (see below). 
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Fig. 2. — WENO numerical solutions of eq.Q with 7 = and b = 0.03. The source is taken 
to be S' = (f){x) and initial condition J{x, 0) = 0. 

Figure 3 plots the time-evolution of J{x, r) with different parameter b, and a zoom-in 
figure at time r = 10^. The zoom-in figure shows clearly that the integral J J{x, T)dx is 
6-independent. This is consistent with the photon number conservation eq. (fT2l) . It shows 
again that the WENO algorithm is robust. We can also see from the right panel of Figure 
3 that all the curves of In J(x, r) vs. x at r = 10^ and in the range —2 < x < 2 can be 
approximated as a straight line. That is, the width xi of the local Boltzmann distribution 
shown in Figure 3 is equaal to about 2, and it is approximately 6- independent. 

The formation and evolution of the local Boltzmann distribution can be quantitatively 
described by B{t) defined as 

i?(r) = 26/[lnJ(0,r)-lnJ(l,r)], (14) 

where In J(0, r) — In J(l, r) is the slope of the straight line In J(x, r) vs. x for |x| < 1. For 
Gaussian source S{x) = (j){x) = e~^'^ /y/rr, we have -B(O) = 2b, and B{t) approaches 2b at 
large r. Figure 4 presents the numerical relation of B{t) vs. r. The slopes [log J(0,r) — 
log J(l,r)] at r = 10^ are, respectively, 0.0601 for b = 0.030, 0.0303 for b = 0.015, 0.0159 
for b = 0.0079, and 0.0051 for b = 0.0025. That is, within the frequency range |x| < xi and 
xi = 1, the relative deviation of the slope din J{x, T)/dx from 2b is no larger than 2%. Thus, 
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r = 10^ can be considered as the time scale of forming a local Boltzmann distribution within 
< xi = 1. For small width xi < 1, this time scale is lower as ~ 10^. Therefore, the time 
scale of the onset of W-F coupling with the width xi equal to about Doppler broadening is 
10^-10^ 



2400 r 




Fig. 3. — WENO numerical solutions of eq.Q with 7 = and b = 0.03 (solid), 0.015 
(dashed), 0.0079 (dot-dot-dashed) and 0.0025 (dot-dashed). The source is taken to be S' = 
(f){x) and initial condition J[x, 0) = 0. The right panel is a zoom-in of the dashed square of 
the left panel. 

We can relate the width xi with the mean number of scattering, Nc, needed to form the 
local Boltzmann distribution. Although the redistribution function eq.(jl]) is 6-dependent, 
the probability of x photons undergoing a resonant scattering per unit time is (j){x), which 
is 6-independent. Thus, at a given time r, the mean number Nc of resonant scattering of 
photons within |x| < Xi approximately is 

1 r' 

Nc^T— / ^(x)dx. (15) 
xi Jo 

Eq.(15) gives the "finite but large number of scattering" for realizing a local Boltzmann dis- 
tribution within X < xi{t) (Wouthuysen, 1952). Therefore, the approximate 6-independence 
of xi (Figures 3 and 4) would imply the 6-independence of Nc. 
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Fig. 4.— B{t) vs. r for solutions of Figure 3 with b = 0.03 (solid), 0.015 (dashed), 0.0079 
(dot-dot-dashed) and 0.0025 (dot-dashed). 



Considering an expanding background, i.e. 7 7^ 0, we solve eq.(9) by the WENO 
algorithm. Figure 5 plots solutions with the same source S = (f){x) and parameter b = 0.03 as 
in Figure 2, but with 7 = 10~^ and 10"^. Similar to Figure 2, a local Boltzmann distribution 
has formed when r > 10'^ for both 7 = 10^'^ and 10^^. The section of the spectrum near 
X = becomes r-independent when r > 10^ for 7 = 10~^, and r > 10^ for 7 = 10~^. We 
call this r-independence to be saturation of the profile J(x, r) around resonant frequency. 
In saturated state, the number of photons redshifted from u > uq to the local Boltzmann 
distribution area ~ z/q due to Hubble expansion is equal to the number of photons leaving 
from uq to the red wing. Therefore, we see from Figure 5 that once J reaches the saturation 
state, the boundary on the red wing [x < 0) of J is moving to left (red). On the other hand, 
the boundary on the blue wing (x > 0) is almost time- independent. 

Unlike in Figure 2, the width xi does not always increase with time. For 7 = 10^'^ the 
width stops to increase when r > 10^, and for 7 = 10~^, it is stopped at r > 10^. One 
can find the mean scattering number with the similar way as eq. lITS!) . When 7 7^ 0, the 
time duration of photons staying in the frequency space from x to x + Ax roughly is Ax/j 
[eq. llTT]) ]. On the other hand, the mean probability of |a;| < xi photons being scattered in a 



unit r is — / (j){x)dx. The larger the x/, the less the probability. Thus, all photons within 



5. Wouthuysen-Field coupling in an expanding background 



5.1. 



Width of the local Boltzmann distribution 




- 10 - 



|x| < xi averagely undergo Nc scattering give by 



1 

Nc ^ — (f)(x)dx. 
7 Jo 



(16) 



From Figure 5, the maximum width for 7 = 10^^ is estimated as x/ = 1.9, corresponding to 
Nc ~ 0.5 X 10^. While for 7 = 10~^, maximum width is xi = 2.8, and A^^^ ~ 5.0 x 10^. Once 
the width xi stops to increse, all quantities in eq. (fT6!) . 7, xi and are r-independent. 

Thus, Nc should also be r-independent. 
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Fig. 5. — WENO numerical solutions of eq.Q with b = 0.03 and 7 = 10~^ (two top panels) 
and 7 = 10^^ (two bottom panels). The source is taken to be S" = 0(x) and initial condition 
J{x, 0) = 0. The right panels are zoom-in of the dashed square of the corresponded left 
panel. 

The r-independence of the width xi is also shown in Figure 6, in which we still use 
7 = 10~^, and J{x,0) = initially. However, the source is taken to be S{x) = (j){x — 10). 
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That is, the source photons have frequency a; = 10, or z/ = z/q + lOAuo- The resonant 
scattering ai x = {u = uq) will occur when these photons have redshifted from u to u = uq, 
which takes time of about r = IO/7 ~ 10^. Figure 6 shows that the whole distribution of 
J{x, t) dramatically evolves with time, but the width of the local Boltzmann distribution 
around x = is kept to be x/ ~ 2 from r = 10^ to 4 x 10^. We also find from our numerical 
calculations that when r > 4 x 10^, the intensity of the photon flux around x = keeps 
constant, or it is in saturated state. 




Fig. 6.— WENO numerical solutions of eq.® with 7 = lO'^ and b = 0.0079 (left), 0.015 
(middle), 0.03 (right). The source is taken to be = (j){x — 10) and initial condition 
J{x, 0) = 0. The bottom panels are the zoom-in of the dashed square of the top panels. 

Similar to Figure 3, Figure 6 also shows that the width of the local Boltzmann distri- 
bution is approximately 6-independent. From eq. fll6p . one can also expect that the width 
will be smaller for larger 7. A local Boltzmann distribution can form only if 7"^ is large 
enough. This property is shown with Figure 7, in which we use the same photon source 
S{x) and parameter b as in Figure 6, but we take larger 7. Figure 7 presents the results of 
7 = 1 and 10~^. We see from Figure 7 that in the case of 7 = 1, there is no local Boltzmann 
distribution at any time. The resonant scattering leads only to a valley around x = 0. It 
is because the strongest scattering is at x ~ 0, which moves photons with frequency x ~ 
to X 7^ 0. However, the redshift lets photon quickly leaving from x ~ 0. They are not 
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undergoing enough number of scattering to form a local statistical equilibrium distribution. 
For 7 = 10~^, it seems to show a local Boltzmann distribution, but its width is very small 
at all time. 




Fig. 7. — WENO numerical solutions of eq.Q with 7 = 1 (left) and 7 = 10^^ (right). 
Parameter b = 0.0079 (dot-dot-dashed), 0.015 (dashed), 0.03 (solid). The source is S* = 
0(x — 10) and initial condition J{x, 0) = 0. The bottom panels are the zoom-in of the 
dashed square of the top panels. 



5.2. Photon source and W-F coupling 

The r- and 6-independencies of the shape and width of the local Boltzmann distribution 
yield an important conclusion that for given parameters 7 and b, the formation and evolution 
of the local Boltzmann distribution is independent of the photon sources S{x,t). This is 
because eq.dH]) is linear of J. Any source S{x, r) can always be considered as a superposition 
of many monochromatic sources around frequency x = Xi, or S{x, r) = J2i Si4>{x — Xi,T), Si 
is the intensity of photon source 0(x — Xi,T) with frequency Xi. J{x, r) can be decomposed 
into J(x, r) = J2i Jii^, t), where Ji{x, r) is the solution of eq.([9]) with the source i. Thus, if 
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the formation of the local Boltzmann distribution around x = is independent of r and b, the 
superposition J{x, r) = Ji{x, r) should also show the same local Boltzmann distribution 
around a; = 0. Although the overall amplitude does depend on the source, the shape around 
resonant frequency does not. 

As an example, Figure 8 presents a solution with the same parameters as in Figure 6, 
but the source is with continuous spectrum given by 



S{x, ■ 



(lO/xf, 10<a;<15, 
otherwise 



(17) 



Photons with frequency a; = 10 will arrive earlier at a; = with higher intensity, while those 
with frequency x = 15 will arrive at x = later with lower intensity. The flux J{x, r) of 
Figure 8 has very different shape from Figure 6, while the local Boltzmann distribution at 
—2 < X < 2 of Figure 8 is exactly the same as that in Figure 6. Therefore, the W-F coupling 
is always working regardless the original spectrum of the redshifted photons. 





Fig. 8.— WENO numerical solutions of eq.© with 7 = 10"^ and b = 0.0079 (dot-dot- 
dashed), 0.015 (dashed), 0.03 (solid). The source is given by eq. (fT7|) . The right panel is a 
zoom-in of the dashed square of left panel. 



5.3. Intensity 

From Figures in §5.1 and 5.2, we see that the intensity of photon flux J at the local 
Boltzmann distribution is strongly dependent on 7 and r. Figures 5 and 6 show that at 
early time J is smaller than its saturation state. For the solution of Figure 6, the flux in the 
frequency range of the local Boltzmann distribution is saturated at about r = 4 x 10^ with 
saturated flux J ~ 10^, while the intensity at r = 10^ is significantly lower than that of the 
saturated state. 
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Figure 9 is the same as Figure 6, but taking 7 = 10"^ and 10~^. The lower panels of 
Figure 9 shows once again that the time-evolution of the intensity is about 6- independent. 
Figure 9 shows also that for 7 = 10~^, the photon flux approaches the saturated state with 
intensity of J ~ 10^ at r = 10^. While for 7 = 10""^, the saturated state has not yet been 
approached even when intensity J ~ 10^, and r = 10^. Generally, the smaller the 7, the 
larger the saturated intensity and the longer the r needed to approach its saturated state. 




Fig. 9. — WENO numerical solutions of eq.Q with 7 = 10 ^ (left) and 10 ^ (right) with 
source S{x) = 0(x - 10) and b = 0.0079 (dot-dot-dashed), 0.015 (dashed), 0.03 (solid). The 
bottom panels are the zoom in of the dashed square of top panels. 

6. Conclusions and discussions 
6.1. Summary 

The onset of the W-F coupling, or the formation of local Boltzmann distribution is 
similar to the process of approaching a statistically thermal equilibrium state via collisions 
or scattering. The particle distribution in the statistical equilibrium is independent of time, 
initial distribution, and the details of collision. The equilibrium distribution is maintained 



15 



only by the enough coUisions among particles. A local Boltzmann distribution is formed 
once the number of resonant scattering is large enough. Like other statistically thermal 
equilibrium, the features of the local Boltzmann distribution are independent of time, photon 
source, initial photon distribution, and etc. 

In an expanding universe, photons are moving in the frequency space with "speed" given 
by the redshift. The formation of the local Boltzmann distribution depends on the compe- 
tition between the resonant scattering and the redshift. If photons have undergone enough 
scattering during their path through the frequency space around the resonant frequency, a 
local Boltzmann distribution will be formed. Otherwise local statistical equilibrium cannot 
be approached. 

In our work, we use the Gaussian profile eq.([3]), but not the Voigt profile. Our major 
results on the formation of local Boltzmann distribution will not be affected by using the 
Voigt profile if the width xi is not larger than the Doppler thermal broadening. We found 
that the solution of eq. (9) with rest background is not affected by Voigt profile if the ratio a 
between the natural and Doppler broadening is equal to 10"'^, which corresponds to T ~ 650 
K. For large a, or small T, the Doppler thermal broadening is small. In this case, the width 
xi, in which the color temperature Tc of Lya photons is locked to the kinetic temperature T 
of hydrogen atoms, would also be small. 



A basic problem for the 21 cm signal from the ionized and heated region around first 
stars is the conditions on which one can estimate the 21 cm emission and absorption with 
the W-F coupling, which forces the internal (spin-) degree of freedom to be determined by 
the thermal motion of the atoms. In this case, the relative occupation of the two hyperfine- 
structure components of the ground state depends only upon the shape of the spectrum near 
the Lya frequency. Therefore, we need a local Boltzmann distribution of Lya photon with 
frequency width equal to or larger than {i^f — i^ol > 1^21 = 1420 MHz, or 



Thus, from Figure 7 and eq.(lO) one can conclude that in regions with /hi < 10~^, where 
7 > 10~^, the W-F coupling will not work. Although electron-hydrogen and proton- hydrogen 
collisions can be important; the 21 cm signal will be incredibly small. On the other hand, 
in the primarily neutral IGM, W-F coupling is very efficient. 



6.2. Applications to the 21 cm problem 




(18) 
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Prom eqs.(7) and (10), we have 




1/2 



10 



) 



3/2 




1/2 



t = 0.26 X l{fh 



,51,-1 



7T yrs. 



(19) 



l + z 



We know that the saturation, or time independent solution around resonant frequency, can 
be used only when the time r is larger than the time of photon moving over a frequency 
space from —Xi to Xi. From Figures 6 and 9, we have r7 equal to about 10 at the saturation. 
Therefore, eq.(19) yields that the time independent solution is available only if the time scale 
of the evolution of the ionized sphere of first stars is larger than about 3 x 10^ yrs. This 
gives a constraint on the 21 cm emission region, as such regions are very narrow, the time 
scale of the evolution being comparable to 10^ yr, or even less (Liu et al. 2007). 

Actually we may not need a saturation state. What we used for estimating the 21 cm 
signals is the frequency distribution J(a;, r) to show a Boltzmann-like shape in the central 
part \x\ ~ 0, i.e. the onset of the W-F couphng. The time-scale of the W-F couphng onset 
is equal to about 10^ yrs for neutral hydrogen density of the concordance ACDM model 
[eq.(7)]. This time scale is much less than that of the evolution of 21 cm region of first stars. 
It seems to indicate that we can safely use the W-F couphng in the 21 cm estimation of first 
stars. 

However, we should mention the effect of the photon intensity. The coupling coefficient 
between Lya photons and spin temperature is proportional to the intensity J (e.g. Furlan- 
etto, Oh, & Briggs 2006). The W-F coupling would generally be suppressed due to the fact 
that the flux at the resonant frequency J{x = 0) is always less than the flux at other frequen- 
cies. In a saturated state this suppression is small (Figure 6). However, before approaching 
the saturated state, the intensity J generally is significantly less than its saturated value. 
That is, although the local Boltzmann distribution is formed at the time of the order of 
T ~ 10^, the intensity at that time would still be low, and the W-F couphng is not enough 
to produce the deviation of Tg from Tcmb- Therefore, it may not be always reasonable to 
assume that the deviation of Tg from Tcmb is mainly due to the W-F coupling if first stars 
or their emission/ absorption regions evolved with time scale equal to or less than Myrs. 

The WENO algorithm revealed the time evolution of photons undergoing resonant 
scattering, whose information is generally lost in the asymptotic solutions, or the time- 
independent solution. Although time-independent solutions provide useful guidance, they 
do not show the conditions for the efficiency of the W-F coupling at different times. The 
asymptotic sohition probably is never reached for short life-time objects. It would be im- 
possible to correctly estimate observable 21 cm signal from ionized and heated halos of first 
stars without a correct understanding of the time evolution of the W-F coupling. 
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The computational domain in the case of static background is x e [—15, 15]. The initial 
condition is J{x, 0) = and the boundary condition is J{x, r) = at the boundaries. In 

the case of expanding background (7 7^ 0), the computational domain is bigger, depending 
on the value of the Sobolcv parameter 7. The domain, denoted as {xieft,Xright), is chosen 
such that J{xieft,T) ~ and J{xright,T) ~ 0. For example, the domain is taken to be 
X e [-100, 15] for the case of 7 = 10"^. 

The computational domain {xieft, x right) is discretized into a uniform mesh as following. 



where N = Ni + Nr and Ax = {xright — xieft)/N, is the mesh size. We also denote Jf = 
J{xi,T"'), the approximate solution values at Xi and the n*'* time step, i.e. = nAt, At 
being the numerical time step. 

A. 2. The WENO algorithm: approximations to spatial derivatives 

To calculate we use the fifth order WENO method (Jiang and Shu, 1996). That is. 



flux in the fifth order WENO approximation because the wind direction is fixed (negative). 
First, we denote 



A. Numerical algorithm 



A.l. Computational domain and computational mesh 



Xi — iAx, 



i^-Nr--,Nr 




hi = J(xi,T"), 



(A2) 



where n is fixed. The numerical flux from the WENO procedure is obtained by 



where h^j^l/2 ^-re the three third order fluxes on three different stencils given by 



(A3) 
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- (2) 1 , 5 1 

"-1+1/2 = 3 * 6 ~ 6 

and the nonlinear weights cUm are given by, 

where e is a parameter to avoid the denominator to become zero and is taken as e = 10^^. 
The hnear weights 7; are given by 

3 3 1 

^^ = 10' ^^=5' ^^=10' ^^'^ 

and the smoothness indicators jSi are given by, 

13 1 
Pi = — (/ii-i-2/i, + /ii+i)2 + -(/ii_i-4/i, + 3/i,+i)2, 

13 1 
(32 = —{hi - 2hi+i + hi+2f + -{hi - hi+2f , 

13 1 



A. 3. Numerical integration: an 0{N) algorithm 
We need to numerically integrate j Tl{x, x')J{x' ,t)dx' , denoted as 

/(x) = i y e2^^'+^'erfc(max(|x + b\, \x' + b\))J{x', t)dx', (A6) 

with 7l{x,x') as in eq.(jlj). To evaluate = I{xm), Vm = —Ni,---,Nr, we apply the 
rectangular rule, which is spectrally accurate for smooth functions vanishing at boundaries, 

1 P-^right 

= - / erfc(max(|x„ + b\, \x' + b\))e^^'''~^'' J{x', t)dx' (A7) 

Nr 

~ ^ eiic{ui8JL{\xm + bl\xi + b\))e^^''^^^^J{xi,t). (A8) 



Notice that this summation algorithm is very costly as it takes 0{N) operations per m, 
therefore the total procedure has 0{N'^) operations overall. We use a grouping technique. 
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described below, so that the overall computational cost can be reduced to 0{N), without 
changing mathematically the algorithm and its accuracy. 

The proposed scheme with order N computational effort is the following. Let Ni, — 
floor {-^) and A^26 = fl(^or{-^). The integration algorithm is designed for two cases: m > 
—Nb and m < —Nb- 

In the case of m > — A^;, or equivalent ly Xm + b >0: 

Im = -Ax{ erfc(|a;, + 6|)e2''-*+^V(x,,t) 

i=-Ni 

m 

+erfc(|a;„ + 6|) J] e^'^^+^" J{xi,t) 

i=-m-N2b 

Nr 

+ ^M\^i + b\)e^'"'^'"j{x,,t)) (A9) 

j=m+l 

= ^Ax(7i,„ + erfc(|x^ + 6|)72,m + /3,m) (AlO) 

1. Evaluate h-N^, h-N^ and h-Ni, respectively as 

h,-N, = J2 eM\xi + b\)e'''^^+'"j{xi,t), (All) 

i=-Ni 

h,.Mb = E e'''''^''j{^r.t). (A12) 
i=Nb-N2b 

Nr 

h^^Mb = erfc(|x, + 6|)e2''-^+''V(x,,i), (A13) 

i=-Nb+l 

which leads to 0{N) cost. 

2. Do m = -Nb + 1, Nr 

Evaluate /i,^, /2,m, -^3,m respectively by 

h,m = /i,m-i - erfc(|x_^_Ar,, + 6|)e'''^— ^26+''V(x_^_Ar,„ (A14) 

/2,m = h,m-i + e^''^-+''V(x„, t) + e^^^— t) (A15) 
/3,m = /3,m-i - erfc(|x„ + 6|)e2^^-+^V(x„, (A16) 

ENDDO 
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To be consistent with the indexes, if Ni — N^b < Nr then, we will set h^rn — 0, for 
m — Ni—N2b, Nj.. The algorithm leads to 0{1) cost per m, therefore 0{N) computation 
overall. 

In the case of m < — iV^, or equivalent ly x^, + 6 < 0: 

^ m— 1 

Im = -Ax{J2 ^M\x^ + b\)e'''''^'''J{x^,t) 

i=-Ni 

-m-N2b-l 

+erfc(|a;^ + 6|) J] e^^^^^^" J{xi,t) 

i=m 

Nr 

+ J2 evk{\xi + b\)e''-^+'"j{xi,t)) (A17) 

i=— m— iV2(, 

= ^Ax{Ii,m + erfc{\x^ + b\)l2,m + h,m) (Al8) 

1. Evaluate /i -atj,-!, h,-Nt,-i and h-Nb-i as 

7i,_^,_i = J2 eM\xi + b\)e^'-^+'"j{xi,t), (A19) 

i=-Ni 

h,-N,-i = Yl e''^'^''j{xi,t), (A20) 

i=-Nb~l 

l3,-N,-i = J2 erfc(|a;, + 6|)e2^^^+^V(a:,,t), (A21) 

i=Ni,+l-N2b 

which leads to 0{N) cost. 

2. Do m = -TVfe - 2, -AT; 

Evaluate /i,^, /2,m, -^3,m respectively by 

= - erfcdx^ + b\)e^''^"^+'^J(x^, t) (A22) 

/2,m = hm+i + e''^-+''j{x,n, t) + e^'^—^^-^+'" J{x.m-N,b-u t) (A23) 
h,m = h,m+i - eTic{\x^m-N,b-i + fe|)e''"— ^26-i+''V(x_„_Ar,,-l, t) (A24) 
ENDDO 

To be consistent with the indexes, if Ni — A26 > Nr, we will set /s,^ = 0, for m — 
—N2b — Nr, —Ni. Again, the algorithm leads to 0{1) cost per m, therefore 0{N) 
computation overall. 
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A. 4. Time evolution 

To evolve in time, we use the third-order TVD Runge Kutta time discretization (Shu 
& Osher, 1988). For systems of ODEs Ut = L{u), the third order Runge-Kutta method is 

= ^ti" + i(7x« + ArL(M«,T" + AT)), 

U^+^ = Jxi"+^(li(2)+ATL(li(2),T" + ^AT)). 
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